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ABSTRACT 

This paper presents the properties of a family of scale-free triaxial haloes. We adduce 
arguments to suggest that the velocity ellipsoids of such models are aligned in conical 
coordinates. We provide an algorithm to find the set of conically aligned velocity 
second moments that support a given density against the gravity field of the halo. 
The case of the logarithmic ellipsoidal model - the simplest triaxial generalisation 
of the familiar isothermal sphere - is examined in detail. The velocity dispersions 
required to hold up the self-consistent model are analytic. The velocity distribution of 
the dark matter can be approximated as a triaxial Gaussian with semiaxes equal to 
the velocity dispersions. 

There are roughly twenty experiments worldwide that are searching for evidence of 
scarce interactions between weakly-interacting massive-particle dark matter (WIMPs) 
and detector nuclei. The annual modulation signal, caused by the Earth's rotation 
around the Sun, is a crucial discriminant between WIMP events and the background. 
The greatest rate is in June, the least in December. We compute the differential 
detection rate for energy deposited by the rare WIMP-nucleus interactions in our 
logarithmic ellipsoidal halo models. Triaxiality and velocity anisotropy change the 
total rate by up to ~ 40%, and have a substantial effect on the amplitude of the annual 
modulation signal. The overall rate is greatest, but the amplitude of the modulation 
is weakest, in our radially anisotropic halo models. Even the sign of the signal can 
be changed. Restricting attention to low energy events, the models predict that the 
maximum rate occurs in December, and not in June. 

Key words: dark matter - galaxies: haloes - galaxies: kinematics and dynamics - 
galaxies: structure - celestial mechanics, stellar dynamics 



The construction of velocity distributions for triaxial haloes 
is a hard problem. Jeans' theorem guarantees that the dis- 
tribution function depends only on the isolating integrals 
of motion. Generally, motion in a triaxial potential admits 
only one exact integral of motion, the energy or the Jacobi 
constant for the case of figure rotation. As Schwarzschild 
(1981) has articulated, it seems that self-consistent triaxial 
equilibrium configurations can exist only if the potential has 
additionally two effective, non-classical integrals. To date, 
the most general and successful method for building tri- 
axial models has been the numerical combination of orbit 
densities, often called Schwarzschild's (1979, 1982) method. 
Recent applications and extensions of the method to triax- 
ial modelling are given by Merritt & Fridman (1996), Zhao 
(1996) and Hafner et al. (2000). 

Given the difficulty of the task, it is often the case that 
only the lowest order velocity moments of the distribution 
function are calculated via the Jeans equations (sometimes 
called the stellar hydrodynamical equations). This is because 
the moments are easier to obtain and are related directly to 



observable properties. Many solutions of the Jeans equa- 
tions have been derived for spherical and axisymmetric sys- 
tems (e.g., Binney & Mamon 1982; Bacon 1985; Fillmore 
1986; Amendt & Cuddeford 1991). Even though it is not 
always evident whether such solutions correspond to dy- 
namical models with positive definite distribution functions, 
these studies have been useful. For example, starting config- 
urations for Ai'-body experiments can be generated by as- 
suming that the velocity distributions are Gaussians with 
semi-axes as prescribed by the Jeans solutions. Of course, 
this is only an approximation, but it is an excellent one and 
in routine use (e.g., Barnes 1994). 

The simplest model of all for a dark halo is the isother- 
mal sphere. This is commonly used to estimate the rates in 
both microlensing and non-baryonic dark matter detection 
experiments (e.g., Paczyhski 1986; Jungman, Kamionkowski 
& Griest 1996; Lewin & Smith 1996). There is no reason 
whatsoever why dark matter haloes should be spherical, and 
there is ample evidence from external galaxies that haloes 
are generically fiattened (e.g., Sackett et al. 1994; Oiling 
1995, 1996). In this paper, we present a number of proper- 
ties of one of the simplest triaxial halo models - the logarith- 
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mic ellipsoidal potential (e.g., Binney & Tremaine 1987; de 
Zeeuw & Pfenniger 1988; Miralda-Escudc & Schwarzschild 
1989). This is the natural triaxial generalisation of the sin- 
gular isothermal sphere. Self-consistent distributions of ve- 
locities for this model are not known, but here we provide 
simple and analytic solutions of the Jeans equations. These 
solutions are aligned in conical coordinates and provide the 
velocity second moments required to support the triaxial 
dark halo against gravity. 

One possibility is that galaxy haloes arc composed of 
weakly-interacting massive-particle (WIMP) dark matter. 
Direct detection experiments measure the energy deposited 
by the rare interactions between the WIMP and the detec- 
tor nucleus. There are now about 20 direct detection exper- 
iments running or in preparation worldwide (e.g., articles 
in Klapdor-Kleingrothaus & Ramachers 1997 and Spooncr 
1997). One of the main difficulties in the experiments is how 
to distinguish between the recoil events caused by WIMPs 
and those caused by radioactivity in the surroundings and 
by cosmic rays. One suggestion is that the annual modula- 
tion in the WIMP signal, caused by the motion of the Earth 
around the Sun, may be a powerful diagnostic (Preese, Prie- 
man & Gould 1985). In fact, one of the experimental groups 
(DAMA) has very recently claimed detection of this modu- 
lation (Bcrnabci ct al. 1999a, b), though the validity of this 
claim is still a matter of fierce debate (e.g., Gerbier et al. 
1997, 1999; Abusaiadi et al. 2000). This paper investigates 
the dependence of this differential rate and the annual mod- 
ulation on the triaxial shape and velocity distribution of the 
dark matter halo. 

The paper is arranged as follows. Section 2 introduces 
the scale-free halo models under scrutiny, while Section 3 
derives and solves the Jeans equations for the velocity sec- 
ond moments under the assumption of conical alignment. 
A number of reasons are given in support of our assump- 
tion as to the orientation of the velocity dispersion tensor. 
The properties of the logarithmic ellipsoidal halo models - 
the intrinsic and projected shapes and velocities - arc found 
in Section 4. Finally, Section 5 considers the application to 
direct detection experiments in some detail. 

2 TRIAXIAL HALO MODELS 

One of the most widely-used axisymmetric halo models is 
the logarithmic scale-free potential investigated by Toomre 
in the 1970s (see Toomre 1982) and subsequently studied by 
others (e.g., Richstone 1980; Evans 1993). Its equipotential 
surfaces are spheroids. The analogous scale-free triaxial halo 
models have equipotentials that are triaxial ellipsoids with 
axis ratios p and q (e.g., Binney 1981; dc Zeeuw & Pfenniger 
1988). We shall call them the logarithmic ellipsoidal models. 
In Cartesian coordinates, the potential-density pair is given 
by: 



^ix,y,z) = ^vl\n{x^ + y'^p ^ + z'^q ^) 



and 



vl Ax^ + By'^p-'^ + Cz'^q- 



with 

C = l+p-^-, 



B = l + q- 



(2.1) 



(2.2) 



(2.3) 




constant fx 
constant v 



Figure 1. The coordinate curves (ju, J^) on the surface of the 
sphere. If the radius of the sphere is r, then the set (r, ix, u) pro- 
vides a triply-orthogonal set of coordinate surfaces. The filled dot 
indicates the point n = v = p^; it lies in the (a;, 2;)-pIane, at 
X = rVr, 2 = rVl -T, with T = (1 _ q,2) ^.^^ triaxiaUty 

parameter. 



The rotation curve is completely flat with amplitude Vc- 
This potential therefore describes a triaxial halo with a 
flat rotation curve. Without loss of generality, we require 
5^ < P'^ < 1- The spherical limit {p = q = 1) is of course 
the familiar singular isothermal sphere. [If desired, models 
with a finite core-radius Rc can be obtained by adding to 
the term in parentheses in the expression for the potential. 
These models are recognized as the triaxial generalization of 
the axisymmetric power- law galaxies (Evans 1993, 1994)]. 

As will become apparent, there is much advantage in 
working in conical coordinates (r, /u, f) defined as follows. 
The first coordinate r is the distance to the origin, so that 
= x'^ + z'^, where {x,y,z) are the standard Cartesian 
coordinates. The variables /u and v are angles. They are the 
solutions for t of 



y 



1 T — p2 ' T — q 



+ ■ 



0, 



(2.4) 



with p and q constants that satisfy the condition < g < 
p < 1. This is a quadratic equation for r with two real 
roots, which we order so that q^ < v < p^ < fJ- < 1- Some 
properties of these coordinates, including the relation with 
the standard spherical coordinates (r, 9, <f), are given in Ap- 
pendix A (see also Morse & Feschbach 1953). 

Figure 1 illustrates the (/j, v) coordinate curves on the 
sphere. The angular coordinates // and v have a center at 
H = V = p^ , which lies at = and 6 = Of, with 



= T = 



1-/ 
l-g2 



(2.5) 



where T is the triaxiality parameter (e.g., Franx, lUingworth 

& dc Zeeuw 1991). While each choice of p and q defines a 
conical coordinate system, all choices which lead to the same 
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value of T correspond to the same set of coordinates. 

The logarithmic ellipsoidal models arc the simplest 
members of the more general family of logarithmic scale-free 
models which have potentials of the form: 



(2.6) 



The function g{ij,, v) is arbitrary, and describes the angular 
dependence of the potential. We write the density distribu- 
tion in our models in the general form: 



pir,ii,iy) 



h{n, v) 



(2.7) 



Here h{^,, v) is another arbitrary function that describes 
the angular dependence of the density. When p is the self- 
consistent density, then a = 2 and the functions h{fi, f) and 
g{lJ,,v) are related by the Poisson equation ATvGp = V'^$. 
This case is useful for the description of the dominant, dark 
matter component of haloes. When a > 2, the density law 
describes a tracer population (such as Population II stars or 
globular clusters) that resides in the halo but falls off faster 
than the dark matter. The Laplacian in conical coordinates 
is given in equation (A6). We find: 



h{ii,iy) = 



(p - v)g' 



■|(m- «^)5^ 



+ 9 



'f I, ^^g 



where w{t) = (r— l)(r — p^)(t — q^), and the prime indicates 
the derivative with respect to the argument. We have taken 
4-kG = 1 and Vc = 1. 

The logarithmic ellipsoidal models are just a special 
case with a simple form of g{p, v): 



p2g2 



(2.8 



The associated self-consistent density (2.2) can be written 
in the form (2.7) with a = 2 and with: 



Hi^, ^) = ^-j [2p^'?^(M + 1')- {p^ + +p^q^)pi'] 



(2.9) 



In what follows, we derive Jeans solutions for general scale- 
free densities in logarithmic scale-free potentials in the next 

Section, but return to the self-consistent logarithmic ellip- 
soidal models as our illustrative examples in Sections 4 and 
5. 



3 THE JEANS EQUATIONS 
3.1 Conical Coordinates 

As Binney & Tremaine (1987) point out, there is a funda- 
mental defect with the Jeans equations - namely, that there 
is no equation of state relating the components of the ve- 
locity dispersion to the density. Any method of solving the 
Jeans equations therefore requires assumptions regarding ei- 
ther the shape or the orientation of the velocity ellipsoid. 
And, of course, the usefulness of such solutions depends on 
whether the assumptions are physically motivated. 

We believe that Jeans solutions aligned roughly with 
conical coordinates do correspond to physical distribution 
functions for triaxial scale-free models. There are three rea- 
sons for this. First, the motion in the separable models 



is described naturally in ellipsoidal coordinates, and many 
positive definite distribution functions exist (Statler 1987; 
Hunter & de Zeeuw 1992). At large radii the density in sep- 
arable models becomes scale-free, and the ellipsoidal coor- 
dinates reduce to conical coordinates. Hence, the velocity 
ellipsoid is conically aligned. Second, in the Jeans solutions 
of axisymmetric scale-free models (see de Zeeuw, Evans & 
Schwarzschild 1996; Evans, Hafner & de Zeeuw 1997), ap- 
proximate distribution functions are constructed for solu- 
tions close to alignment in spherical polar coordinates. Con- 
ical coordinates appear to be the most natural generalisation 
of this result to triaxial scale-free models. Finally, Amendt 
& Cuddeford (1991) have argued that Jeans solutions for 
which the kurtosis vanishes can be regarded as physical. For 
triaxial scale-free models, the vanishing of the kurtosis im- 
plies alignment in conical coordinates. These arguments are 
suggestive, rather than rigorous, but they seem worth pur- 
suing to us. 

Let us consider a triaxial halo with a gravitational 

potential $ = ^{r,fi,u) and with a phase-space distribu- 
tion function / = f{r,ij,,i',Vr,v^,Vv). In equilibrium, / 
must satisfy the coUisionless Boltzmann equation (Binney 
& Tremaine 1987, §4.1) which, in conical coordinates, be- 
comes: 



= ru, 



21 

dr 



2 , 2 



df 



dr J dvr 



- VrVf^ 



Of 



- VrV„- 



df 



1 

1 



df 
9f 



.2{ii- 



df 

dp, J dVf, 



v^v„ df 



dv 



+ 



v) dp J dVf, 2{p, - 

d^\df_ VjA 

p) dv ) dv,, 2{i' - 



v) dvv 

. df 



^2{v — p) dv J dvu 2{h' — p) dvf_, 

where Q and R are the metric coefficients given in equa- 
tion (A5). Multiplication by Vr, v^j, and v^, respectively, and 
subsequent integration over velocity space, gives the Jeans 
equations. These relate the mass density p and the elements 
of the stress tensor p{viVj) (with i,j equal to r,p,iy) to the 
forces. We obtain 

dp{vl) , p r^/_.2\ 2\ 2\1 , 1 \dp{VrV^) 



dr r 

p{VrVy.) 



2{vl) - (vl) ~ {vl) 



1 



dp 



-+ 



2{p - v) 



+ 



1 



dp{VrVu) p(VrVv) 



df 



+ 



2{v-p) 



9£ 
'~d^' 



(3.1a) 



dp{vij,v,,} p{Vfj.Vi,) 



dp{vl) 



dv 

P{{vl) - (vl)) 



{v-p) 



dp 



2{p- 



_p_d^ 

rQ dp ' 



dplvrVu) „p , 1 1 
or r rQ 



dp{vvVij) p{vvVij,) 



1 

Vr 



dp{vl) 



dp 

P{{vl) - {vD) 



dv 



(p-v) 

p_a£ 

rR dv 



(3.1b) 



(3.1c) 



2{v - p) 

These three relations between the six stresses p{v^), p(v^), 
p{v1), p{vrVi_i), p{vrVv), and p{vy,Vv) must be satisfied at 
any point in a triaxial halo model with potential $ and 
mass density p. A solution of these equations for given p and 
$ corresponds to a physical equilibrium model only if it is 
associated with a distribution function /(r, /x, v,Vr,Vfi,Vv) > 
0. 

These throe partial differential equations must be sup- 
plied with boundary conditions. The stresses must all vanish 
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at infinity. Furthermore, the factor (/i — v) vanishes at the 



special point fi 



p . In order to avoid singularities in 



the terms on the left-hand side, when /Lt = = p we must 
have (c.f., Evans & Lynden-Bell 1989) 

(vl) = (vl), {Vf^V„) = {VrV^} = {VrV^} = 0. (3.2) 

At this point, the velocity ellipsoid is isotropic in the angu- 
lar direction. However, the radial dispersion (vr) may differ 
from (vft) = (v^) here. 

3.2 The Scale-Free Ansatz 

The potential (2.6) and density (2.7) have the desirable at- 
tribute of scale-freeness. Their properties at radius r' = kr 

follow from those at radius r by a simple magnification, 
and by a rescaling of the time variable t' — kt. For ex- 
ample, p{kr,n,v) = k~" p{r, fijv). We consider distribu- 
tion functions that are also scale- free, i.e., that satisfy 
f{kr,iJ,,v,Vr,v^,Vu) = k'" f{r,ii,i',Vr,Vix,Vv)- The associ- 
ated stresses then have the following form: 



2\ _ F2{fi,iy) 



2^ Fs(p,u) 



p{VrVy,) = 



F4{n, v) 

FiiilJ-, v) 



(3.3) 



where F\, F2, F3, F4, F5 and Fe are functions of n,iy. A 
necessary - but not sufiicient - condition for the stresses 
(3.3) to correspond to a physical equilibrium model is that 
Fi,F2,Fs > 0, since they give the velocity average of the 
non-negative quantities v^, v^, and v^. F4, F5 and Fe may 
be negative, but we always require that the eigenvalues of 
the stress tensor are non-negative. Even then, not all such 
solutions correspond to positive distribution functions. 

Substitution of the forms (3.3) into the Jeans equations 
leads to a system of three partial differential equations for 
the six variables Fi,...,6. We find : 



(a - 2)Fi +F2 + F3 = h + 



1 

+ R 



Q 

dFe 
dv 



dFs 
dfj, 



2(/i 

Fe 



2{v - n) 



dF2 _^ F2-Fi 
d/j, 2{ii — v) 



dFs F3 - 
dv ^ 2{v 



h dg 
2g 9/1 
dFi 



+ Q{a - 3)F5 



+ 



Fa 



dv {v — ji) 



M) 



h_dg_ 

2g dv 



+ R{a - 3)F6 



dF\ ^ F4 
dfj, ifJ- — I 



(3.4a) 



(3.4b) 



(3.4c) 



We shall find it convenient to refer to the right-hand sides of 
(3.4a), (3.4b) and (3.4c) as Ki{fi, u), K2{iJ.,v) and K-s{fi,h'), 
respectively. We are free to pick three of the six functions 
Fi,...,6 arbitrarily and solve (3.4) for the other three. 

3.3 Solution for Conical Alignment 

Suppose we pick F4, F5 and Fa and solve for Fi, F2 and i*!^. 
This choice is motivated by the fact that it is then easy to 



consider the special case -Fk.s.e = 0, in which the velocity 
ellipsoid is everywhere aligned exactly along the conical co- 
ordinate system. We are particularly interested in this align- 
ment, as we have argued that it is likely to correspond to 
physical distribution functions. 

For the special case a = 2 (the dark halo case) , a simple 
approach can be used. In this ceise Fi drops out of equation 
(3.4a). We use the remaining relation to rewrite (3.4b, c) as 
follows: 

dF2 _^ F2 _h r 1 _ Idg- 
2 Ifj, — v g djjL. 
ft r 1 1 dg'\ ^ C3 
2\_v — ti g dv . 



dFa 

dv V — fj, 



+ 



ji — V 
F3 



+ 



C2 



2(m- 



.V — n 

where the functions C2 (/i, v] 



(3.5) 



2it,-v)' 
and 03(11, v) are defined as 



C2 



C3 



F4 + {V- //) 



dFj 
dv 



Q 



2{n - v)QF5 + 



R 



dFs^ 
dp, 
dFe 
dv 



Fs 



2{p - v) 

Fe 



+ 



F4 + {n- V 



dF\ 
dp. 



dFr, 



+ 



2{v - p) 

F5 



dp 2(p — v) 



2{v - p)RFe + 



R 



dFe 



+ 



Fe 



dv 2{y — p) 



(3.6a) 



(3.6b) 



Equations (3.5a) and (3.5b) can be integrated separately, by 
using p — V and v — p, respectively, as integrating factors. 
The result is 



F2 = 



2{p- 



2{v~p) 



G2{p,v) + C2{p,v) + D{v) 
G3{p,v) + C3{p,v) + D{p) 



(3.7) 



where, for the moment, the functions D{v) and D{p) are 
arbitrary, and we have introduced the functions 



(v - m) dg{m,v) i 
dm ^ 



G2 = dmh(m,v)\l+ , ^ 
J ^ 9{m,v) 



(3.8) 



and 



dmc2(m, v), 



-I 



dnc3{p,n), (3.9) 



so that G2 and G3 depend only on the chosen density and 
potential. Equation (3.7) gives the general solution for F2 
and -F3 when a = 2, subject to the requirement on the func- 
tions D{p) and D{v) imposed by the boundary conditions 
(3.2), and by equation (3.4a). This latter requirement can 
be written as: 

D{p) - D{v) = 2{p - v)h G3 - G2 + Cs - C2. (3.10) 

The boundary condition (3.2) then shows that D{p^) = 

D{p^), so that we can in fact consider D and D the same 
function A(t), with A{p) = b{p) and A{iy) = D{v). When 
F4, Fs and Fe are chosen such that the right-hand side of 
this equation can be written as a difference A{p) — A{v), 
then the entire Jeans solution is specified: Fi is also arbi- 
trary (but non- negative), and F2 and F3 are given by (3.7). 
Otherwise, one can specify only two of the three functions 
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F4, Fs, and Fe,, and attempt to solve equation (3.10) for the 
third. For ^4,5,6 = both C2 and C3 vanish, and the solutions 
without cross terms in the stress tensor have C2 = C3 = 0. 

When a ^ 2 (the tracer population case), it is best 
to differentiate (3.4b) with respect to f, differentiate (3.4c) 
with respect to fj, and subtract, to obtain the following equa- 
tion for the difference F*(/u, u) = F2{tJ,, v) — Fslfj,, u) : 



d^F* 



+ 



(---) 



^dK,_u^ (3.11) 



djj, 



The right-hand side of (3.11) contains only known functions. 
This equation for F* (n, v) is identical to an equation solved 
by Evans & Lynden-Bell (1989) by the method of Green's 
functions. If the initial data F*{p^,v) and F*{pL,q^) is pro- 
vided, then F*{ii,v) is given everywhere as a double inte- 
gral over the Green's function (Evans & Lynden-Bell 1989, 
section 4). Once we have the difference F*{fjb,i/), then it is 
straightforward to solve for the individual components: 



F2(m, 



F2{j>\iy)+ / dm / . 

J n 2(m - v) 



(3.12) 



With F2{i-i,v) and Fo^{^,v) now available, F\[^,v) is ob- 
tained by rc-arranging (3.4a). Therefore, the boundary con- 
ditions required are: (i) the difference F*{jj,, v) = J2(a*, — 
Fsin, v) on the surface fJ- = and on v = q^, together with 
(a) the value of either F2{n, v) or F3{n, v) on either ji = 
or on V = . 



4 THE LOGARITHMIC ELLIPSOIDAL HALO 

4.1 Intrinsic and Projected Shapes 

In this section, we return to the logarithmic ellipsoidal mod- 
els - partly for their astrophysical importance, partly for 
their simplicity. The density law (2.2) is sensibly positive 
whenever 



q^ +p'^q^ -p^ > 0. 



(4.1) 



The physically allowed domain is illustrated as the hatched 
region in Figure 2. The oblate models have = 1, while 
prolate models have = 1^- Although the equipotential 
surfaces are ellipsoidal, the density figures deviate from a 
pure ellipsoidal shape. Cross-sections of the intrinsic shape 
with the three principal pianos for the model with p — 0.9 
and q = 0.8 are shown in Figure 3. The Eixis ratios of the 
density distribution are 



02^ ^ /(pV-fp'-g') 

ai J q^ -\- p'^ — p^cp- ' 

g'(pV + g'-p') 



V ffli / g2 -h — 



(4.2) 



p^q^ 



For our example with p = 0.9 and q — 0.8, the density con- 
tours have semiaxes in the ratio 1 : 0.774 : 0.489. The flat- 
tening in the density is of course greater than the flattening 
in the potential. 

In order to calculate the projected shape, we choose new 
coordinates {x" ,y" , z") with the z"-axis along the line of 
sight, and the a;"-axis in the {x, j/)-plane (see e.g., de Zeeuw 
& Franx 1989). The z-axis of the triaxial model projects 
onto the y" axis, and x" and y" arc Cartesian coordinates 
in the plane of the sky. The various transformations needed 




Figure 2. The logarithmic cUipsoidal models have density distri- 
butions that arc everywhere positive in a hatched region of the 
(p^, g^)-plane. Oblate models have = I, while prolate models 
have = q^- 

to calculate the projected surface density are summarised in 
Appendix B. The required integration reduces to an integral 
given in Appendix C of Evans & de Zeeuw (1994). The pro- 
jected surface density E(a;",y") is analytic for all viewing 
angles {'&,ip), namely: 



[x ,y )=npq 



x"^ + 



{CIX"^ - C2X"y" + C3t/"^)3/2 



where 

ci = sin^ ifi -\-p^ cos^ ifi, 

C2 = 2(1 — p^) sin 1^ cos ip cos 

2 2 2 2 2 2 2 

C3 = cos ip cos -& +p sin tp cos -& + q sin 1 



(4.3) 



(4.4) 



When the direction of observation lies in one of the three 
principal planes (i9 = n/2,ip — 0,ip = 7r/2), the coefficient 
C2 vanishes, and the principal axes of the projected surface 
density lie along the x"- and ^/''-directions. However, for all 
other viewing directions C2 7^ (for triaxial models with 
p 7^ 1), the minor axis of the projection is misaligned from 
the projected short axis of the model, which by definition 
falls along the 2/"-axis. Since our models are scale-free in 
projection, the misalignment is independent of radius. Fig- 
ure 4 gives contour plots for the projected surface density of 
the model with p = 0.9, q = 0.8 for three different viewing 
directions. 

Let us introduce polar coordinates {R', @) in the sky 
plane, defined by the relations x" = —R' sin®, y" = 
R' COS®. So, the position angle © is defined with respect 
to the j/"-axis and is measured in the counter-clockwise di- 
rection. The projected surface density is E oc S{Q)/R'. The 
resulting formula can be simplified further by introducing 
the two parameters 

p. =|(ci -I-C3) + |V(c3 -ciy+4, 
V =|(ci -I- C3) - |\/(C3 - ci)2 -I- c|. 

These are the conical coordinates (/i, of the direction of 
observation defined by the angles (t?, ip) (c.f. eq. [5.4] of de 
Zeeuw & Franx (1989)). We obtain: 



S(i?',e) 



"Tipq 



R'\jx + v-{jl-v) cos(2e - 2e,)]3/2 ' 



(4.5) 
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Figure 3. Contour plots of the intrinsic density in the principal planes for the logarithmic ellipsoidal model with p = 0.9 and q = 0.8. 




Figure 4. Projected surface density for the logajrithmic ellipsoidal model with p = 0.9 and q = 0.8, and three different directions of 
observation, namely a) i9 = 0° and = 0°, b) i? = 30° and = 45°, and c) = 90° and v = 90°. 



where we have defined the misalignment angle G* by 



90° (see Figure 4). 



tan 20* 



C2 



C3 - Cl 



(4.6) 



If Ot satisfies equation (4.6), then so docs 0, + n/2. If 
we choose 0, to be the root for which (cs — ci)cos26, + 
C2sin2Gt < 0, then G, is the position angle of the major 
axis, as measured counter-clockwise from the j/"-axis. Our 
expression (4.6) for the misalignment angle is identical to 
that found for the classical ellipsoids, in which the density 
is stratified on similar concentric ellipsoids with semi-axes 
1 : p : q (c.f. Stark 1977). Since the potentials (2.6) are strati- 
fied on such ellipsoids, the major axis of the projected poten- 
tial (and hence the projected surface density) must lie along 
a position angle 0*. The axis ratio b'/a of the isophotes, 
defined by the condition E(fe', G* — -1) = E(a', G*), is given 
by 



a' \p.J 



3/2 



(4.7) 



If p = 0.9 and q = 0.8, the projected surface density has a 
shape roughly like E3 if i? = lys = 0° and like E5 if i9 = <p = 



4.2 Velocity Second Moments 

We now carry through the algorithm of Section 3 and find 
a simple set of velocity dispersions that support the loga- 
rithmic ellipsoidal model. Since a = 2, Fi is arbitrary, while 
F4 — F5 = Fe all vanish for conical alignment. Upon substi- 
tution of the expressions for g{fj,, v) and /i(/i, v), the integrals 
in (3.7) can be carried out and we obtain: 



F2{^l,v) = 



2{^l - v) 



[p +q +p q )- 



(4.8a) 



2{u - p) 



/ 2 , 2 , 22 
(p +q +p q 



2pV)l 



(4.8b) 
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The condition F2 + F3 = h can be satisfied by the following 
choice: 



(4.9) 



with A either the function D or the function D, r respec- 
tively either v or /x, and ^ is a constant. The requirement 

that F2 and F3 arc finite at the point ^ = v — means that 
the term in square brackets must vanish there. This means 
that we must choose ^ = 0. The expressions for F2 and ^3 
can then be simplified to 



(4.10) 



Despite the symmetry of F2 and F3 with respect to /x and 
v, the velocity ellipsoid is anisotropic everywhere, except at 
the point ^ = v = p^. Here F2 = F3, while Fi remains 
arbitrary. 

To gain insight into a physically reasonable choice for 
Fi, let us consider the spherical limit {p = q = 1). The 
logarithmic ellipsoidal model becomes the well-known sin- 
gular isothermal sphere. The angular stresses corresponding 
to (4.10) become 



(vl) = {vl) = i, 



(4.11) 



This is recognised as the Jeans solution generated by the 
constant anisotropy distribution functions for the isothermal 
sphere (see eq. [5.10] in Evans (1994)). The full solutions for 
the spherically aligned stresses are 



{vl) = 



1 



2 + 7' 



{vl) = {vl) = 



(4.12) 



where 7 is a constant anisotropy parameter. In the spherical 
limit, 7 = gives an isotropic solution, 7 — » — 2 is the radial 
orbit model, while 7 — > 00 is the circular orbit model. This 
suggests a possible choice for F\ as 



Fi = 



2-^7 (2-|-7)p2g2 



(4.13) 



Any choice for Fi satisfies the Jeans equations. This choice 
(4.13) has the advantage that the spherical limit {p^ = = 
1) certainly corresponds to a positive definite distribution 
function for any choice of anisotropy parameter 7 > — 2. The 
inclusion of the term g{n,i') in (4.13) additionally ensures 
that the axisymmetric limit {p^ = 1) with 7 = has the 
same stresses as the two integral distribution function (see 
eq. (2.3) of Evans (1993)). It is clearly heartening that the 
checkable limits do correspond to physical models. Even so, 
the tria^xial models probably only have physical distribution 
functions for some range in 7 near to the isotropic value. 
Strongly anisotropic models are known to be afflicted by 
instabilities, like the radial orbit instability (e.g.. Palmer 
1994). Henceforth, we only consider models in which the 
inequality in any two semiaxes of the velocity dispersion 
tensor is no greater than 3:1. 

Given our solutions for F\,F2 and Fz, the intrinsic ve- 
locity second moments of the logarithmic ellipsoidal model 



are found through (3.3) as; 



{vl) 



(2 + 7)j32g2 [2(^ + y)p2qZ _ (p2 _^ g2 + p2g2)^^] ' 

_ 1 (3At + u)p\^ - ip' + + p'q^)^lv 



=i 

^ 2 2(/i + t^)p2f/2 _ (p2 + g2 _|_ p2,j2^f^i, ' 



(4.14) 



K) 



1 + ■iz/jp-fi- - [p- + q- + p-q-)lii-' 



2 2(/x + v)p'^q'^ — (p2 + q2 _|_ p2g2^^^ • 



The relationship between conicals (/i, v) and familiar Carte- 
sian coordinates is given in eqs. (Al) and (A2). This fur- 
nishes an analytic and realistic sot of stresses that support 
the logarithmic ellipsoidal model - a prototype for a tricix- 
ial dark matter halo. This is one of the main results of the 
paper. 

The logarithmic ellipsoidal model, much like the isother- 
mal sphere itself, does not obey the virial theorem (e.g., Bin- 
ney & Tremaine 1987, Gerhard 1991). In fact, the velocity 
dispersions (4.14) have a curious property - the total kinetic 
energy is not fixed, and the radial velocity dispersion may be 
changed independently of the angular velocity dispersions. 
The same property holds for the isothermal sphere itself, 
as is evident from (4.12). It is a surprising feature, as we 
expect increases in the radial velocity dispersion to be bal- 
anced by decreases in the angular dispersions. In fact, both 
isothermal spheres and ellipsoids make complete sense when 
regarded as the inner parts of finite mass models. For ex- 
ample, the truncated, flat rotation curve model (Wilkinson 
& Evans 1999) has a density that varies like within a 
scalelength a so that: 



p(r) 



r2(r2 -Ka2)3/2' 



$(r) = — log 



\/r2~P"a2' + a 



,(4.15) 



where, as for the logarithmic ellipsoidal model, we are using 
units with Vc = AttG = 1. The radial velocity dispersion is 
(see eq. (11) of Wilkinson & Evans (1999)) 



1+2 (r^ + a^) 



2\3/2 



Jo 



^"'^ (a) a3 cosh^ 



3inh^+T , 
cosh^ 6 



(4.16) 



In the inner parts, the rotation curve is flat with unit am- 
plitude and the model looks like the isothermal sphere. In 
the limit r << a, the crucial point is that the velocity dis- 
persions (4.16) do indeed reduce to those of the isothermal 
sphere (4.12), as a careful Taylor expansion demonstrates. 
This is the case for all 7 > —2. The kinetic energy of the 
inner, isothermal parts of these models can be increased - 
as it can be balanced by a loss in the kinetic energy of the 
outer parts. 



5 PARTICLE DARK MATTER DETECTION 
RATES 

5.1 Preliminaries 

The dark matter in galaxy haloes may be composed of 
weakly-interacting massive particles (WIMPs), which cou- 
ple to ordinary matter only through electroweak-scale in- 
teractions. Particle physics has provided an extensive list of 
candidates (e.g., Kolb & Turner 1989), of which the lightest 
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Figure 5. The distribution of speeds and the difFciential rate on the major axis {a,b) and the intermediate axis (c,d) of the logarithmic 
cUipsoidal model with p = 0.9 and q = 0.8. In each panel, results arc given for a radially anisotropic [7 = —1.78] solution (dotted line) 
and tangentially anisotropic [7 = 16] solution (dashed line), as well as a comparison Maxwellian (full line). Panel (b) assumes that the 
WIMPs are scattering off '^^Ge nuclei, while panel (d) assumes that they are scattering of -^^Na nuclei. The computations are carried out 
for the date of June 2nd when the total rate is at a peak. 



stable neutral supersymmetric particle (generally the neu- 
tralino) is one of the current favourites (e.g., Jungman et 
al. 1996). One promising way of confirming this hypothe- 
sis involves direct detection experiments. Broadly speaking, 
the experiments work by measuring the recoil energy of a 
nucleus in a low background laboratory detector which has 
undergone a collision with a WIMP. The aim is to measure 
the number of events per day per kilogram of detector ma- 
terial as a function of the recoil energy Q. Although this 
deposited energy is minute and the WIMP-nucleus interac- 
tion is very rare, there are a number of such experiments 
in progress around the world. These include the UKDMC 
collaboration operating in Boulby mine (e.g.. Smith et al. 
1996), the DAMA collaboration in the Gran Sasso Laborar 
tory (e.g., Bernabei et al. 1999), which both use NaT scintil- 
lators, and the CDMS experiment located underground at 
Stanford University, which uses cryogenic germanium and 



silicon detectors (e.g., Gaitskell et al. 1997). 

In all these experiments, the detection rate depends on 
the mass and cross-section ao of the WIMP, as well as 
the mass of the target nucleus ttin in the detector. But, it 
also depends on the local dark matter density po and the 
speed distribution fs{v) of WIMPs in the Galactic halo near 
the Earth. Calculations have already been performed us- 
ing Maxwellian velocity distributions for singular and cored 
isothermal spheres, as well as for self-consistent flattened 
halo models (e.g., Frocso et al. 1985; Jungman et al. 1996; 
Kamionkowski & Kinkhabwala 1998; BeUi et al. 1999). One 
of our aims here is to assess the likely uncertainties in 
the detection rates caused by halo triaxiality and velocity 
anisotropy. The formulae for the calculation of rates in di- 
rect detection experiments arc summarised in the review of 
Jungman et al. (1996). We give here only the bare details. 
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Figure 6. The WIMP annual modulation signal is shown for a radially anisotropic [7 = —1.78] model (dotted line), a tangentially 
anisotropic [7 = 16] model (dashed line) and a comparison MaxwcUian (full line). The upper panel gives the variation in the total rate, 
the lower panel the variation in the rate of low energy events (< 10 keV). This calculation assumes 40 GeV WIMPs interacting with 
^^Ge nuclei. The halo model has p = 0.9 and q = 0.8, while the sun is located on the major axis. 



The differential rate for WIMP detection is 



dR 



a"oPo 



dv, 



(4.17) 



where R is the rate, Q is the recoil energy and mr = 
mN?Tix/("^N + w-x) is the reduced mass, F{Q) is the nuclear 
form factor, fs (v) is the probability distribution of WIMP 
speeds relative to the detector, and 



2mr 



1/2 



(4.18) 



The most commonly- used nuclear form factor is (e.g., Ahlen 
et al. 1987; Preese et al. 1988; Jungman et al. 1996) 



where Qo is the nuclear coherence energy 



3h' 



2mN-RN'' ' 

and -Rn is the radius of the target nucleus (in cm) 
= 10'" [0.3 + 0.91(mN/GeV)'/^] . 



(4.19) 



(4.20) 



(4.21) 



The total event rate can be found by integrating over all 
detectable energies 



R 



P —dQ 



(4.22) 



where Et is the threshold energy for the detector. 

As our benchmark model, we take a WIMP with only 
scalar interactions of mass rux = 40 GeV and a cross-section 

(Jo = 4 X 10'""' cm^. Wc consider two kinds of detectors. 
The first is a cryogenic detector made of germanium so that 
TTiN = 68 GeV. The second is a scintillation detector made 
of Nal, for which results are presented for WIMPs scatter- 
ing off sodium nuclei (mN = 22 GeV). Threshold effects are 
neglected, so that Et = 0. The local halo density po is taken 
as 0.3 GeVcm"^. These values axe suggested by Jungman 
et al. (1996) as a standard set. We note that they are, of 
course, subject to very substantial uncertainties. In particu- 
lar, speculations as to WIMP masses range from 10 GeV to a 
few TeV. WIMPs may have both scalar and spin-dependent 
interactions with the nucleus. In fact, direct detection ex- 
periments are best suited for scalar-coupled WIMPs, whilst 
indirect methods (such as searching for WIMP annihilation 
products, like neutrinos, in the Sun) are more powerful for 
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spin-coupled WIMPs (Kamionkowski, Griest, Jungman & 
Sadoulet 1995). Finally, the uncertainty in the local halo 
density is at least a factor of two, and possibly more (e.g., 
Bahcall, Schmidt & Soneira 1983; Turner 1986; Gates, Gyuk 
& Turner 1995) 

5.2 TricLxiality and Anisotropy 

Let us consider two possible locations for the Sun, namely on 

the long axis (a;-axis) and the intermediate axis (y-axis) of 
the triaxial halo. At these locations, the conical coordinates 
are locally equivalent to cylindrical polar coordinates. So, 
the velocity distribution in the Earth's rest frame can be 
approximated as the triaxial Gaussian 



f- 



(27r) 2 (Tua^a, 



exp[- 



20-^ 



.(4.23) 



Here, wo assume that the Earth moves with respect to the 
rest frame of the Galaxy with a velocity 



V(s = 220 



1.05 + 0.07 cos(V^) 

V lyr 



kms 



(4.24) 



where tp corresponds to June 2nd. This gives an annual mod- 
ulation to the WIMP signal, which is of course invaluable in 
distinguishing it from background events caused by radioac- 
tivity and cosmic rays (Freese et al. 1988). Even though 
many of the experiments operate underground to cut down 
the effects from cosmic rays and use high-purity material to 
minimise the effects of radioactivity, it is still true that the 
background rate is a few events per kilogram per day. This is 
larger than that rate expected from WIMP interactions by 
a factor of ~ 10. The WIMP signal is expected to attain a 
maximum in early June and a minimum in early December. 
Recently, the DAMA group has claimed detection of this 
modulation (Bernabei et al. 1999a,b), although their inter- 
pretation has been contested by others (Gerbier et al. 1999; 
Abusaiadi et al. 2000). 

On the major axis, the Jeans solution (4.14) becomes 



ctr' ={vl) 



(2 + 7)(p-2 + g-2-l)' 



2 / 2\ 



2 ^ vl{2q-^-l) 

2(p-2 -I- g-2 - 1) ' 



(4.25) 



2 / 2\ 



^;2(2p-2_l) 



2(p-2 _ ly 

On the intermediate axis, we have 

2 —4 

2,2, VcP 

Or =(vr} = 



2 / 2\ 



(2 + 7)(l + 9-2-p-2)' 

«?(2g-2-p-2) 



2{l+q-^-p~ 



(4.26) 



„ 2 ,2, _ vl{2-p-^) 

-{v^) - 2(l + g-2_p-2)- 



The distribution of WIMP speeds /^(v) follows from (4.23) 
as 

fniv) = 3 / da sin a / df3 

{2-n)2aRa^az Jo Jo 

r (vcosa + Vm)'^ sin^ a sin^ /? t; sin^ q cos^ /3 n 

expl — 07/ _ r- _ / I 



2CT^2 



2aR^ 



2(Tz 



It is now straightforward to calculate the differential rate 
(4.17) for any of our triaxial, anisotropic halo models. 

Figure 5 shows the distribution of speeds and the differ- 
ential rate for the logarithmic ellipsoidal halo with p = 0.9 
and q = 0.8. The Sun is located on the major axis (a;-axis) for 
panels (a) and (b), the intermediate axis (y-axis) for panels 
(c) and (d). The curves are drawn for a radially anisotropic 
and a tangentially anisotropic velocity distribution that sup- 
ports the triaxial figure. Also shown is the standard curve 
obtained by assuming a spherical halo with a Maxwellian 
velocity distribution (c.f., Fig. 22 of Jungman et al. 1996). 
Panel (b) shows the case when the WIMPS interact with ger- 
manium detector nuclei, while panel (d) the case of sodium 
nuclei. The effects of triaxiality and velocity anisotropy can 
cause the total rate to vary by ~ 20% in the case of germa- 
nium and by ~ 40% in the case of sodium. The total rate 
is greatest for the radial anisotropic velocity distributions, 
least for the tangential anisotropic. In the case displayed in 
Figure 5, the radially anisotropic model [7 = —1.78] has 
a broader distribution of speeds. This gives larger incident 
WIMP velocities, and so more events at higher energies. 
The converse is true for the tangentially anisotropic model 
[7 = 16], which generate more events at lower energies. 

More important than the differential rate is the size of 
the annual modulation effect. The upper panel of Figure 6 
shows the variation in the total rate during the course of the 
year for the radially and tangentially anisotropic halo mod- 
els. This is for the case of 40 GeV WIMPs impinging on ger- 
manium nuclei. Also shown is the comparison Maxwellian 
of the isothermal sphere, for which the annual modula- 
tion effect is already small; the peak to peak variation 
over the course of the year is just ~ 0.019 events/kg/day. 
For the radially anisotropic model, this halves to ~ 0.010 
evonts/kg/day, while for the tangentially anisotropic model, 
this increases slightly to ~ 0.024 events/kg/day. Let us re- 
call that the modulation must be detected against a back- 
ground rate that is at least a factor of 10 greater than the 
underlying WIMP event rate. Even if the WIMP velocities 
are described by a Maxwellian, the detection of the annual 
modulation is hard enough. It is disconcerting to find that 
the modulation amplitude can halve on moving to radially 
anisotropic, flattened models. The lower panel of Figure 6 
shows the variation in the rate of low energy events only 
(< 10 keV). The amplitude of the modulation is severely at- 
tenuated. For all the velocity distributions, the peak to peak 
variation is very low, at most 0.008 events/kg/day. Indeed, 
even the sign of the modulation has reversed, and the max- 
inmm number of events now occurs in December, not June 
(c.L, Hasenbalg 1998). 

The DAMA group (Bernabei et al., 1999a,b) have re- 
cently extracted statistical evidence for the annual modula- 
tion signal in the low energy events using nine 9.7 kg Nal 
detectors. It is therefore interesting to carry out the calcu- 
lations for this case as well. Figure 7 shows the annual vari- 
ation in rate for WIMPs scattering off sodium nuclei in Nal 
detectors. Again, we see that the annual modulation in the 
total rate is weakest in the radially anisotropic model, show- 
ing a peak to peak variation of 0.08 cvcnts/kg/day; for the 
tangentially anisotropic model, it is 0.12 events/kg/day. The 
lower panel shows the modulation in the rate of events with 
measured energies less than 10 keV. Nuclear recoils of ener- 
gies in the 1-20 keV range produce little ionisation as they 
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Figure 7. As Figure 6, but the calculation assumes 40 GeV WIMPs scattering off ^"''Na nuclei. The upper panel shows the variation in 
total rate, the lower panel the variation in the rate of low energy events (< 10 keV) assuming a quenching factor of 0.3. 



lose energy; most of the energy goes into phonons. Hence, 
ionisation or scintillation detectors are much less ofBcient 
than cryogenic detectors. Accordingly, we have included a 
"quenching factor" of 0.3 to account for the efhciency of the 
detector (Spooner et al. 1994; Bernabei et al. 1996); a mea- 
sured energy of 10 keV corresponds to a nuclear recoil en- 
ergy of ^ 33 keV. Wc see that the peak to peak variation in 
the low energy events (which is where DAMA have claimed 
evidence of the signal) is sensitive to the anisotropy. The 
lowest variation is just 0.01 events/kg/day for the radially 
anisotropic model; this rises to 0.07 events/kg/day for the 
tangentially anisotropic model. Reversed modulation (i.e., 
the maximum occurring in December) still occurs, but the 
quenching drives it to still lower energy ranges. For example, 
it occurs for events with measured energies less than 5 keV 
in the models with both Maxwellian and radially anisotropic 
velocity distributions. 

6 CONCLUSIONS 

In this paper, we have presented a paradigm for the triaxial 
dark halo. The logarithmic ellipsoidal model is the natu- 
ral generalisation of the isothermal sphere into the triaxial 
domain. Many of its properties are disarmingly simple - in- 
cluding the potential, mass density and projected surface 
density. This paper has provided simple, analytic solutions 



of the Jeans equations for the logarithmic ellipsoidal model, 
under the assumption of conical alignment of the velocity 
ellipsoid. These solutions give the components of the veloc- 
ity dispersion tensor required to hold up the triaxial halo 
against gravity. In this paper, we have adopted the approx- 
imation that the velocity distribution is a triaxial Gaussian 
with somiaxcs equal to the velocity dispersions as specified 
by our solutions of the Jeans equations. While this distribu- 
tion function does not satisfy Jeans theorem (e.g., Binney & 
Tremaine 1987), it does have the right stresses or momentum 
flux to hold the halo up against gravity. The logarithmic el- 
lipsoidal models therefore provide simple prototypes for the 
density and velocity distributions of triaxial haloes. 

Solutions of the Jeans equations contain as a subset 
the models which have physical meaning, i.e., which have 
a non-negative distribution function. Building distribution 
functions in the triaxial case is not easy, since only the en- 
ergy is known to be an exact integral of motion in the gen- 
eral case. Merritt & Pridman (1996), following up an earlier 
suggestion by Schwarzschild (1993), have provided numeri- 
cal evidence that cusped triaxial systems may undergo slow 
evolution towards axisymmetry. This may mean that exact 
distribution functions for such models do not exist - leaving 
the Jeans equations as one of the few possible investigative 
tools for these slowly evolving systems. 

As an application, we have considered direct detection 



12 N.W. Evans, CM. Carollo, P.T. de Zeeuw 



rates of particle dark matter. Weakly-interacting massive 
particle (WIMP) detection experiments typically measure 
the nuclear recoil energy as the WIMPs collide with detector 
nuclei. Predictions for WIMP detection rates are habitually 
carried out using the isothermal sphere to model the dark 
matter distribution. This paper has provided calculations for 
more realistic triaxial dark haloes with anisotropic velocity 
distributions. Standard isothermal sphere calculations give 
estimates for the total rate that are good to within ~ 20% for 
Ge detectors and to within ~ 40% for Nal detectors. Both 
the shape of the differential rate distribution and the size 
of the annual modulation effect are sensitive to the velocity 
distributions. In pessimistic cases, the modulation ampli- 
tude roughly halves on moving to radially anisotropic and 
flattened models, rendering detection of this characteristic 
WIMP signature even harder. In optimistic cases, tangen- 
tially anisotropic models can give slightly higher peak to 
peak changes than standard isothermal spheres, although 
the total rate is smaller. Assuming a WIMP mass of 40 GeV 
and a cross-section ctq = 4 x 10~^^ cm^ for Ge nuclei, then 
the peak to peak variation in the total rate over the year 
may then be as low as ~ 0.010 events/kg/day or as high 
as ~ 0.024 events/kg/day, with the maximum occurring in 
June and the minimum in December. 

If only the low energy events are considered, then the 
amplitude of the modulation is smaller still and even the sign 
of the correlation can reverse. The greatest number of low 
energy events can occur in December, not June. The energy 
at which this crossover happens depends on the anisotropy 
and triaxiality of the halo, as well as the detector mate- 
rial. In practice, direct detection experiments monitor the 
number of events within energy intervals chosen by the ex- 
perimenter. It would be interesting to calculate the optimum 
energy ranges for detection of both the modulation in the 
high energy events and the reversed modulation in the low 
energy events. In such difficult experiments, where the back- 
ground is typically ten times or more as strong as the signal, 
this effect should be exploited to provide a convincing and 
unambiguous signature of actual detection. 
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APPENDIX A: CONICAL COORDINATES 

Here, wc collect some properties of the conical coordinates 
{r,fi,v) defined in Section 2.1. The two roots /j, and v of 
equation (2.4) are 



or, inverting the relations, by: 
xvx + yvy + zvz 



Vr 



2R 



in - 



+ (/i - 1)(m - q'jyvy + (/i - 1)(a* - p^)zv.] ,(^8) 



+ q^)yvy + (z/ - l)(z/ - 



ti,V=\{kx±^2), (Al) 

with 

r^fci = (j/^ + z^) + p^{x^ + z^) + q^ix^ + y'^), 

/fc2 = lip' - q')x' + (1 - q')y' + (1 - p^)z^f (A2) 

It follows that 11 + 1' = ki and A/xv = kf — are each simple 
rational functions of x'^ , y^, and The inverse transforma- 
tions are: 



^2_ r'{n~q'Ki^-q') 



(A3) 



so that each point (r, v) corresponds to eight points 

(±x, ±j/, iz). The standard spherical coordinates (r,6,4>) 
are related to the conical coordinates by 



tan d) = 



(I_q2)(p2_q2)' 

(M^P^)(p^-^)(l-g^) 

{1 - H){1 - i^){p^ - q^) ■ 



(A4) 



The metric of the conical coordinate system is given by 
ds^ = dx'^ + dy"^ + dz^ = dr'^ + r'^iQ^dp,'^ + R^du^), with 



2 _ I' - 



4w(/u)' 4w(z/) 
u;(r) =(r - l)(r - p^)(r - g^). 



(A5) 



The clement of area dA on the sphere of radius r and the 
Laplacian are given by: 



dA=r'QRd^,d^^ ' (f^--)'^l^d. 



4:y/-w{p)y/w{u) 



^2 1 d 2d 1 d R d 1 d Q d 
V = ——r— + ^r---^ — — — + 



(A6) 



J.2 Qj. J.2QJI dnQdn r'^QR dvRdv' 



The relations between the conical and cartesian velocity 
components can be written as: 



^--^L«'^+2Q(^_l) + 2Ji(!.-l)J' 

Vj, = - Vr- + 

r 

= - \vr + 



r 2Q(^-p2) 2i?(i/-p2) 



(A7) 



+ 



2Q(m-52) 2i?(z.-g2)J' 



APPENDIX B: THE PROJECTED MOMENTS 
OF TRIAXIAL MODELS 

Here, we derive expressions for the projected moments of a 
triaxial model. The projected zeroth moment is the surface 
density, calculated in the main body of the paper. Formulae 
for the projected first and second velocity moments are also 
given. 

We choose new coordinates (a;", y" , z") with the z"-axis 
along the line of sight, and the a;"-axis in the (a::, j/)-plane 
(e.g., de Zeeuw & Franx 1989). This means that the z-axis 
projects onto the j/"-axis. We define the direction of obser- 
vation by the viewing angles {'&,<f), and write r = {x,y,z). 
Then the coordinate transformation is 



Rsr, 



where 



R3 



— sm if 
- cos y cos i) 
cos ifi sin 1? 



R— 1 // 



cos ^ 

- sin ifi cos •& sin i? 
sin (fi sin & cos 6 , 



(Bl) 



(B2) 



( — sin — cos <p cos i? cos <^ sin i? ^ 
cos — sin cos •& sin ip sin •& 
sini9 cosi? , 

The projected surface density is S where: 
E(x",j,")= 1 d^"p(R3'r"). 



(B3) 



The major and minor axis of the projected surface density E 
of a triaxial galaxy generally do not line up with the x"- and 
j/"-axes, but lie at a position angle G* =^ 0. The expression 
for 0, is given in cq. (4.6). Let (x ,y') be Cartesian coordi- 
nates on the plane of the sky aligned with the principal axes 
of the projected density. Then; 

r" = Rar', r' = R^^r", (B4) 

where 

/ cos 0, — sin 0, \ 
R2 = I sin0* COS0* , 











(B5) 



/ cosG. sinG. 0' 
Ro^ = - sine, cosG, 











For completeness, let us also outline the calculation of 
the projected velocity moments. By differentiation with re- 
spect to time, we find for the velocity components: 



r = Rsr, 



1-,— 1.// 
r = R .O r . 



(B6) 
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At a fixed position, we can calculate the intrinsic mean mo- 
tion by integrating over all velocities weighted by the nor- 
malised velocity distribution function. Subsequent integra- 
tion along the line of sight, weighted with the density dis- 
tribution, then gives the observed mean motions: 

1 f°° II II 1 Z""" 

(Mx"> =^ / dz" p{x"), inyii) = - dz" p{y"}, 

(vios) =(Pz") = ^ / dz' p{z'), 

J — OO 

The projected second moments follow in a similar way; for 
example, the diagonal components of the second-rank tensor 
are 



1 

{l^l"x") =^ / dz" p{x"'^), 

■J —oc 



{lh"y")=-^j ^z"p{y"^), (B8) 

J —OO 

(vL) = {l4"z>') =^ / dz" p{z'^). 

and similarly for the mixed components. The products of 
the velocity components in the integrand can be found by 
use of equation (B6). 
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